Maximal Mean Field Solutions in the Random 
Field Ising Model: the Pattern of the 
Symmetry Breaking 

Marco GU AGNELLI^), Enzo MARINARI^' 5 ) 

and 

Giorgio PARISI^ 

(a): Dipartimento di Fisica and Infn, 
Universita di Roma Tor Vergata, 
Viale della Ricerca Scientifica, 00173 Roma, Italy 

(6): Department of Physics and NPAC, 
Syracuse University 
Syracuse, NY 13244, USA 

(a): Dipartimento di Fisica and Infn, 

Universita di Roma La Sapienza, 
P. Aldo Moro 2, 00185 Roma, Italy 

guagnelli@roma2 . infn. it marinariOromal . infn. it parisiOromal . infn. it 



March 20, 1993 



Abstract 

In this note we study the mean field equations for the 3d Random 
Field Ising Model. We discuss the phase diagram of the model, and 
we address the problem of finding if such equations admit more than 
one solution. We find two different critical values of (3: one where the 
magnetization takes a non-zero expectation value, and one where we 
start to have more than one solution to the mean field equation. We 
find that, inside a given solution, there are no divergent correlation 
lengths. 
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1 Introduction 



The Random Field Ising Model (RFIM) (see for example refs. [1, 2, 3]) 
is waiting for pieces of new understanding and further clarifications of the 
relevant physical mechanisms. 

Let us start by sketching the theoretical situation. For a certain time 
it was hoped that dimensional reduction could be the appropriate method 
to compute the critical behavior of a ferromagnet in presence of a random 
magnetic field. It was proven in [4] that in perturbation theory the sum of 
the most divergent diagrams close to the phase transition for a random field 
model in dimension D coincides with that of a ferromagnetic theory, without 
random field, in the reduced dimension d = D — 2. The terms that are 
neglected are less singular than the leading ones by a factor £~ 2 , £ being as 
usual the correlation length. This result suggests that all the exponents of the 
random field system coincide with those of the corresponding ferromagnetic 
system in D — 2 dimensions. 

Clearly this result cannot be correct. Simple physical arguments (con- 
firmed by a rigorous analysis [5]) lead to the conclusion that the lower critical 
dimension is 2, not 3, as implied by dimensional analysis. The deep reason 
for this failure can be found following the non-perturbative analysis of ref. 
[3, 6]. Let us summarize the main results. 

We assume that the system is described by the following Hamiltonian 
density, which is a functional of the order parameter <f>(x): 



where the random field h(x) is a Gaussian uncorrected white noise with 
variance g 5(x — y), and g parametrizes the strength of the random field. 

The stationary points of H can be found by solving the corresponding 
mean field equations 



When these equations admit only one solution, as it happens for suffi- 
ciently large temperature, it is natural to introduce the correlations functions 




(1) 



+ V'{(j>) = h(x) . 



(2) 



C{x) 



0(x)0(O) , 
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(3) 



where by the long bar we denote the thermal average over all the realizations 
of the random magnetic fields. 

These two correlation functions are the mean field approximation to 
< 4>(x) >< 0(0) > and < <j)(x)<j)(0) > c respectively. Then one finds that C(x) 
is proportional to the same correlation function of the pure system in dimen- 
sions d = D — 2. The functions G(x) and C(x) are related one to the other. 
In Fourier space one finds that 

C(k) = gfd»Mjp±^ 

G(k) = Jd^ piu) ^J-j . (4) 

The function G(k) is the same as for the pure system in dimension d = 
D — 2 (dimensional reduction works in configuration space with the function 
C, and in momentum space with the function G). 

When 2 admits more than one solution to compute expectation values we 
must assign a weight to each solution. This makes life more complicated. If 
we label by a different solutions, and by w a the relative weight we can write 



C{x) = ]>> a 0(:r)<KO) 



G \ X ) = l^ W c 



5h(0) 



By using different prescriptions for the weights w a we can obtain different 
results. This is especially true if the number of different solutions of the mean 
field equations increases with the volume. 

Dimensional reduction can still hold, but with a crazy choice of the 
weights: 
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_ sign detj-A + V"^)] 

w a — 7 , \V) 

where Z w is such that J2a w a = 1- Here the Morse theorem states that 
Z w = 1. 

A physically motivated choice would be: 

sign det[-A + V"{<t> a )\ exp(-(3H [0 a ]) 

W* = 7} , (7) 

where a runs over all the solutions of the mean field equations, minima, 
maxima and saddle points all together. The strange looking unusual factor 
sign det[— A + V"(<f> a )] is needed to keep the continuity of Z w when new 
solutions appear. 

It could be argued that the energy of the minima is so smaller than the 
energy of the saddle points and maxima, that we can simply write 

_ exp(-f3H[(f)a]) 

w a — 7 i V->) 



and keep the sum restricted only to the minima. In the rest of this paper we 
follow this second strategy. 

It is possible that this modified mean field theory gives the correct results 
(as it is implicit in the work of ref. [7]) and that the failure of dimensional 
reduction is simply related to the existence of many solutions with different 
energy ([3, 6, 8]). 

Our aim is here to investigate numerically this improved mean field ap- 
proximation to make its predictions explicit and eventually to compare them 
with Monte Carlo simulations. We have been motivated to start this inves- 
tigation by an interesting paper [9], in which it was suggested that replica 
symmetry is already broken at the point ferromagnetic phase transition. For 
results obtained both in the mean field framework and with a Monte Carlo 
and a T = optimization approach, see refs. [10, 11, 7]. 

In this note we limit ourselves to the study of two particular solutions of 
the mean fields equations, which we call 0+ and 0_. They are such that for 
any solution <p a (and for any x) the relation (f)-(x) < (f) a ( x ) < (fi+( x ) holds. 
The existence of two solutions with this property (in the high temperature 
phase they coincide) follows from convexity arguments [12]. We call them 
maximal mean field solutions. 
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2 Lattice Mean Field Equations 



We consider the Random Field Ising Model (RFIM) with Ising type (Z 2 ) 
variables defined on a 3c? simple cubic lattice. We study the solutions of its 
mean field equations. 

With i we denote the triplet of integers (x,y,z), which characterize the 
lattice sites. We will consider configurations of the random field {hi = 
Oi H}, where the quenched random variables 6i can take the values ±1 with 
probability |, and we have chosen the absolute value of the field, 7i, to be 
1.5. Such choice for Tt was meant to allow the critical temperature T c to 
have a non negligible shift from T c in the pure model, and simultaneously 
not to be large enough to allow the transition to become first order ^. 

In the mean field approximation one introduces local magnetization vari- 
ables TOj, which play the same role of <f>(x) in the continuum formalism. The 
total free energy is written as a function of the local magnetization, and the 
condition for the free energy being stationary is the usual mean field equation 

mi = tanh(/3 (Vrrii + hi)) , (9) 

where with Dwij we define the lattice sum over the 6 first neighbor variables. 

If this equation admits only one solution there is no ambiguity. If, on 
the contrary, there are many solutions, one has to weight (according to the 
previous discussion) different solutions with a weight proportional to the 
exponential of minus the free energy (multiplied by /3). 

Our ideal goal is to look for all solutions of this equation, which correspond 
to local minima of the free energy, but this is an awful task when the number 
of solutions is very large, as it happens at low T. Here we have just looked 
for the solutions with higher, positive and negative, magnetization, (ra+ and 
m_) using a simple iterative scheme. We have started the iterative procedure 
used to solve eq. (9) from the two initial conditions to, = m s and to* = — m s . 
Although a completely safe procedure would start from m s — 1, it is more 
convenient (and it does not change the results) to take a value for m s slightly 
smaller than one. The appropriate value of m s depends on the temperature; 
in our simulations we have taken m s = .6. 

In the high T regime both these runs converge to the same (unique) 
solution. In a broken phase they will tend to different solutions with average 
magnetization of opposite signs. This procedure should be good enough to 
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localize the temperature T below which the solution of (9) is not unique, and 
to give relevant quantitative hints about the structure of the phase transition. 

We label the solutions of the mean field equations, in a given realization of 
the magnetic field, by the index a; given the pattern of our search a is limited 
to take only one or two values. For each realization of the magnetic field the 
index a belongs to the set A (which can be, in our simulation, constituted of 
1 or 2 solutions). The average over different field samples (which we denote 
by a bar: we denote the average over different solutions by (•)) is done by 
having A running from 1 to Na- 

In each solution a (characterized by the V = L 3 values of the local 
magnetization rrii) we compute the relevant observables. We define the total 
magnetization density 



V 

and the sum of the squared local variables 



™ Q = iE<, (io) 



We define the energy density 



E° = -±Z(±m?1»n? + him?) , (12) 
the entropy density 

^4D^¥^) + ^M^», (13) 

and the total free energy as 

F a = V{f3E a - S a ) . (14) 
The weight w a associated to each solution a is given by 

exp(-/?F«) 

Wa = = • (15) 
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3 Numerical Results for Local Quantities 



Here we present numerical results for system of size up to 48 in a range of (3 
that goes from 1.1 to 1.5 (we will always give f3 in units of the critical /3 at 
zero random field, i.e., |). We have analyzed 600 random field samples for 
the 12 3 lattice, 400 for the 24 3 lattice, 200 for the 36 3 lattice and 30 for the 
48 3 lattice. 

In this section we will discuss the behavior of local quantities (i.e., those 
objects that can easily be constructed from the magnetization), while in 
the next section we will concentrate our attention on the response functions, 
which must be computed by inverting the lattice equivalent of (—A + V"(4>)), 
a highly non-local operation. 

A very interesting quantity is 

w* = Y,<- (is) 

a 

This quantity is different from 1 when the mean field equations admit more 
than one solution: roughly speaking W~ 2 is the average number of relevant 
solutions. We display the results for W 2 as function of f3 in fig. 1. We see 
that W 2 becomes sizably different from 1 only at /3 greater than 1.2. We see 
a change in regime at this beta, which we denote by 

Another quantity that is interesting to measure is the maximal magne- 
tization m 2 M , defined as max Q (m a ) 2 . In fig. 2 we show the f3 dependence of 
m 2 M for different lattice sizes. We see a transition from an asymptotic zero 
value of m 2 M to a non zero value around (3 = 1.35. The transition becomes 
sharper by increasing the size of the lattice. We see a change in regime also 
at this new value of (3, which we denote fa- 

A more detailed understanding can be obtained by considering the cor- 
relation functions of the local magnetization. At this end we define, for each 
solution a, the magnetization on a 2-plane as 

M:(X) = J2m a (X,y,z) . (17) 

My(X) and M"(A) are defined in an analogous way. We define the zero 
(bi-) momentum magnetization-magnetization correlation function for the so- 
lution a as 
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C a (\)= E M-(X 1 )M^(X 2 ) . (18) 

[l=X,y,Z] Ai,A2 such that Ai — A2 1 =A 

The total correlation function at distance A, averaged over N A samples, 
is defined as 

A a 

and we denote by C C (A) its connected part. 

At first order in perturbation theory in the strength of the random field 
[3, 6, 4] C C (A) has (as we have discussed before) a double pole in Fourier 
space. It has also been shown that this form retains its validity at all orders in 
perturbation theory, and should be exact in the region where supersymmetric 
considerations hold. In x space that leads to 

C C (A) ~ A(l + 7^)e"«^ + B , (20) 

that defines the correlation length 

In fig. 3 we plot the inverse correlation length as a function of (3. We 
have used a global fit to C(A) (which has determined A and B, by as- 
suming a functional dependence that takes in account the periodic boundary 
conditions). In all cases we have computed the statistical errors by using a 
standard jack-knife procedure. We have also computed A dependent corre- 
lation length estimators. By averaging them in the plateau region we have 
obtained another estimate of which turns out to be completely compat- 
ible with the one coming from the global fits. The fits turn out to be of very 
good quality, confirming the approximate validity of the form (20). 

The correlation length of fig. 3 has quite a broad maximum close to /3 — 
1.35. close to its peak increases when going from L = 12 to L = 24, but 
for larger lattices it remains constant. 

In fig. 4 we plot the maximum value of the correlation length, as a 
function of -jr, to stress the saturation that occurs for large L. The asymptotic 
is of order 4.5. It is rather consistent that the correlation length becomes 
independent from the size for sizes 3 to 4 times larger than the correlation 
length. 
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In fig. 5 we plot the coefficient B (i.e., the constant asymptotic value 
of the correlation function) computed from the fit to C(A) as a function of 
(3. For large volumes B should become identical with m 2 (which, in our 
analysis, turns out to be very similar to m 2 M ), but its finite size corrections 
are smaller, especially in the high temperature region, where B and /77^^ £1X6 
asymptotically zero. B seems to take a non-zero expectation value starting 
from (3 2 — 1.35. This method gives a very good estimate of the value of the 
critical temperature where m 2 M becomes sizably different from zero. 

In the region where m 2 M is zero all different solutions of the mean field 
equations should become locally equal in the infinite volume limit, or more 
precisely their absolute difference should be in average go to zero with the 
volume. 

It is natural to ask if these values of (3 do correspond in the thermody- 
namical limit to real phase transitions. The quantity W 2 becomes different 
from zero as soon as there exist a realization of the magnetic field that ad- 
mits two solutions. An explicit computation shows that if h(i) = (— \Y +y+z 
one finds two solutions when (3 > (3q — 1.015 1 . Simple minded arguments 
(which generalize the original Griffiths theorem for random diluted magnetic 
systems) suggest that the free energy is C°° but not analytic at (3q 

For (3 > (3g there exist realizations of the magnetic field for which at least 
two solutions exist. These field configurations are special, and their measure 
is small. We expect therefore that W 2 is different from 1, mathematically 
speaking, for (3 > (3q, but it becomes sizably different from 1 only at f3 > f3\. 
Similar arguments can be done for (3 2 . The non vanishing of m 2 M for (3 > (3q 
is a pathology that arises from our choice of considering only the maximal 
solution. If we consider the physically relevant quantity, i.e., 

m 2 = $> a (m a ) 2 , (21) 

a 

it should become different from zero only at values of (3 much higher than 
(3q. The fact that the correlation length remains finite and somewhat small 
near f3 2 may be taken as an indication that the true ferromagnetic transition 
at which m 2 becomes different from zero is at higher values of (3. 

The situation would be clarified if we could compute the full expression 
for C(x), summing over all the solutions, but we have left this task for a 

1 An approximate formula valid for small TL is (3q = 1 + + 0(H 4 ). 
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future work. 



4 Numerical Results for the Response Func- 
tions 

To compute the correlation functions in the mean field approach we must 
use the fluctuation-dissipation theorem. We are therefore lead to consider 
the susceptibility function Xi,j, which is equal to the derivative of the mag- 
netization rrii with respect to the field hj (for sake of typographical clearness 
in the following we will omit the solution label a). If there is a single stable 
state we have to perturb the unique solution of equation (9). We get in this 
way the equation 

Xi j = 0(1 - m?)(Z>Xij + Stj) . (22) 

This is a linear sparse equation that can be solved by using standard iterative 
techniques. 

The computation of x f° r & U the value of % and j would be extremely time 
consuming, so we compute the Green functions gi = Xi,o by setting j = 
and iterating the relation 

gi = P(l-m*)(V gi + 8 ifi ) . (23) 
We also compute the susceptibility X = y Xi,j by iterating 



Xl = (3(l-mf)(V Xi + l) , 

X = ^Ex,, (24) 
v i=i 

and the overlap susceptibility x q = y Xi,j m i m j from 



x\ = Pii-m^VxHrn) , 
1 v 

x q = 77 E Xi m ■ 

i=l 



(25) 
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The name overlap susceptibility arises from the following considerations. 
Let us consider two replicas (a and r) of the same system whose dynamics is 
determined by a Hamiltonian that contains the usual one system contribution 
plus a direct coupling among the two systems. The total Hamiltonian is 

H[a]+H[r}-eJ2^ri ■ (26) 

i 

This construction is common in the study of other disordered systems like 
spin glasses. The quantity \ q coincide with evaluated at e = 0, where q 
is the overlap density, i.e., y J2i o~%Ti- 

In the interesting case in which the mean field equations admit many 
solutions a we follow the simplest procedure of weighting each of these with 
the weight w a (we remind that we are taking in account only the two maximal 
solutions). In this way we are obtaining only one term of the two that form 
the full susceptibility. It is easy to check that the response function 

R(U) = S-E^ (27) 
orij a 

is given by 

R(i, j) = Xij + 0(2 w <* m i m< j ~ w a m i 2 w -r m ] ) • ( 28 ) 

a a 7 

The second term (in brackets), which arises in presence of many solutions, is 
likely to be dominant near the critical point, as will shall see below. It may 
be convenient to call the first the diagonal contribution, and the second one 
the off-diagonal contribution. 

We have computed the diagonal contributions x an d X 9 with the results 
shown in figs. 6 and 7. It is impressive that x has a sharp maximum close to 
Pi, while x q has a peak at much higher beta (slightly above P2) and does not 
show any significant anomaly at f3i . These two peaks are volume independent 
for large volume. The correctness of this result is confirmed by the direct 
analysis of the correlation length corresponding to x, which we show 
m fig. 8. 4 1} does 

never become large in the whole region and for (3 < Pi 
essentially coincides with the correlation length we have discussed in 
the previous section. We find that the supersymmetry prediction of equality 
of the two correlation lengths is correct in the region P < Pi where only one 
solution is present. 
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We have also considered the correlation lengths £^ n ) defined by taking 
the n-th power of the zero bi- momentum correlation functions that give 
and by looking at their decay. They do not present a significant difference 
(once divided by n) from the one obtained for n — 1. In fig. 9 we show the 
correlation length with n = 2 (which has the smallest statistical error), which 
can be compared with the n — 1 result of fig. 8. The two sets of curves are 
very similar. 

Evaluating at least some approximation to the off-diagonal contribution 
to the susceptibility is essential. We have done it by only using our maximal 
solutions. There is a large statistical error. In the low temperature region 
we expect that the off-diagonal contribution is proportional to iV 1//2 , this 
contribution arising from a few exceptional configurations of the magnetic 
field that have two solutions with opposite magnetization with similar weight. 
This event happens with a probability of order 1/iV 1 / 2 ; the corresponding 
off-diagonal susceptibility is of order N, so that the net contribution to the 
susceptibility coming from these exceptional configurations is proportional to 
TV 1 / 2 . In this region the off-diagonal susceptibility is dramatically increasing, 
showing the trend to diverge about j3 ~ 1.30. Anyhow there is no convincing 
argument that implies that the restriction to the maximal solution should 
be a good approximation, apart from very close to fli, where only two stable 
solutions are expected. 

A full computation (including all the solutions) of both the diagonal and 
the non-diagonal contribution to the susceptibility would be extremely inter- 
esting. 

5 Conclusions 

The existence of many solutions to the mean field equations turns out to be 
a crucial phenomenon; inside a single solution (at least of the maximal type) 
one does not see any sign of the presence of a divergent correlation length. 
The critical behavior of the susceptibility and of the correlation length in a 
3d RFIM is dominated by the effects of the presence of many solutions. The 
super symmetric predictions start to fail exactly at the point where one finds 
more than one solution of the mean field equations. It is not surprising that 
dimensional reduction, which completely misses the existence of more than 
one solution, gives unreliable exponents at the critical point. 
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It is reasonable that each solution of the mean field equation does cor- 
respond to a valley for the energy in configuration space 2 . In this case the 
dynamics of Monte Carlo simulations of a real system also at temperature 
slightly above the critical one is likely to be dominated by thermal activated 
tunnelling among different valleys, and we expect it to be a slow process. 
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2 We reserve the word state for solutions (or cluster of solutions) such that their distance 
v Si( m f — m l) 2 docs not vanish in the infinite volume limit. 
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Figure Captions 

1. W 2 as a function of (3. dots for the 12 3 lattice, dashes-dots for the 24 3 
lattice, dashes for the 36 3 lattice and solid line for the 48 3 lattice. 

2. As in fig. 1, but m 2 M . 

3. The inverse correlation length m = as a function of j3 for different 
lattice sizes. 

4. The maximum correlation length as a function of the inverse linear 
size of the system. 

5. The constant coefficient B from the global fit to C(A) as a function of 
P. 

6. As in fig. 1, but x> the diagonal contribution to the susceptibility. 

7. As in fig. 1, but x q i the diagonal contribution to the overlap suscepti- 
bility. 

8. As in fig. 1, but fW. 

9. As in fig. 1, but 4 2) - 
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